Modular Control of Biological Networks

The concept of control is central to understanding and applications of biological network models. Some of their key structural features relate to control functions, through gene regulation, signaling, or metabolic mechanisms, and computational models need to encode these. Applications of models often focus on model-based control, such as in biomedicine or metabolic engineering. This paper presents an approach to model-based control that exploits two common features of biological networks, namely their modular structure and canalizing features of their regulatory mechanisms. The paper focuses on intracellular regulatory networks, represented by Boolean network models. A main result of this paper is that control strategies can be identified by focusing on one module at a time. This paper also presents a criterion based on canalizing features of the regulatory rules to identify modules that do not contribute to network control and can be excluded. For even moderately sized networks, finding global control inputs is computationally very challenging. The modular approach presented here leads to a highly efficient approach to solving this problem. This approach is applied to a published Boolean network model of blood cancer large granular lymphocyte (T-LGL) leukemia to identify a minimal control set that achieves a desired control objective.


Introduction
With the availability of more experimental data and information about the structure of biological networks, computational modeling can capture increasingly complex features of biological networks [1,2].However, the increased size and complexity of dynamic network models also poses challenges in understanding and applying their structure as a tool for model-based control, important for a range of applications [3,4].This is our focus here.To narrow the scope of the problems we address we limit ourselves to intracellular networks represented by Boolean network (BN) models.BNs are widely used in molecular systems biology to capture the coarse-grained dynamics of a variety of regulatory networks [5].They have been shown to provide a good approximation of the dynamics of continuous processes [6].
For the commonly-used modeling framework of ordinary differential equations, there is a well-developed theory of optimal control, which is largely absent for other modeling frameworks, such as Boolean networks or agent-based models, both frequently used in systems biology and biomedicine.Furthermore, control inputs, in many cases, are of a binary nature, such as gene knockouts or the blocking of mechanisms.For BNs, there is no readily available mathematical theory that could be used for control, leaving sampling and simulation.As networks get larger, with hundreds [7] or even thousands of nodes [8], this leaves few computational tools to identify control inputs for achieving preselected objectives, such as moving a network from one phenotype (e.g., cancer) to another (e.g., normal).One approach is to reduce the system in a way that the reduced system maintains relevant dynamical properties such as its attractors [9,10].This allows the control methods to be applied to the reduced system, and the same controls can then be used for the original system.Control targets in Boolean networks have been identified by a variety of approaches: using stable motifs [11,12], feedback vertex sets [13,14], model checking [15,16], and other methods [17,18].A few approaches have used strongly connected components for model analysis by decomposing the wiring diagram or the state space of the network [19,20].However, to our knowledge, none of the approaches developed thus far exploit the modular structure exhibited by many biological systems, in order to identify control strategies by focusing on one module at a time, which is the approach used in this paper.
Modularity refers to the division of the system into separate units, or modules, that each have a specific function [21,22].Modularity is a fundamental property of biological systems that is essential for the evolution of new functions and the development of robustness [23,24].In [25], we developed a mathematical theory of modularity for Boolean network models and showed that one can identify network-level control inputs at the modular level.That is, we obtain global control inputs by identifying them at the local, modular level and assembling them to global control.This enables network control for much larger networks than would otherwise be computationally feasible.It is worth noting that, although the ideas and concepts behind our approach to modularity are intuitive and natural, to our knowledge, there is no well-developed mathematical theory of modular decomposition of Boolean networks and the use of modularity for the purpose of control.In this paper, we further develop this approach into a mathematical theory of biological network control.
We further propose to use another property of biological networks, represented through Boolean network models.Almost all Boolean rules that describe the dynamics of over 120 published, expert-curated biological Boolean network models have the property that they exhibit some degree of canalization [26].A Boolean function is canalizing if it has one or more variables that, when they take on a particular value, they determine the value of the function, irrespective of the values of all the other variables.As an example, any variable in a conjunctive rule (e.g., x ∩ y ∩ z) determines the value of the entire rule, when it takes on the value 0. We derive a criterion for Boolean network models whose Boolean functions are all canalizing, that can be used to exclude certain modules from needing to be considered for the identification of controls.
Our approach to control via modularity is summarized in Figure 1.We decompose the network into its constituent modules, then apply control methods to each module to identify a control target for the entire network.We show that by combining the controls of the modules, we can control the entire network.In the last part of the paper, we present theoretical results that exploit the canalizing properties of the regulatory functions to exclude certain modules from the control search.Finally, we demonstrate our approach by applying it to a published model of the blood cancer large granular lymphocyte (T-LGL) leukemia [27].

Background
We first describe Boolean networks and how to decompose a network into modules.In a BN, each gene is represented by a node that can be in one of two states: ON or OFF.Time is discretized as well, and the state of a gene at the next time step is determined by a Boolean function that takes as input the current states of a subset of the nodes in the BN.The dependence of a gene on the state of another gene can be graphically represented by a directed edge, and the wiring diagram contains all such dependencies.

Boolean Networks
Boolean networks can be seen as discrete dynamical systems.Specifically, consider n variables x1, . . ., xn each of which can take values in F2 := {0, 1}, where F2 is the field with two elements, 0 and 1, where arithmetic is performed modulo 2.Then, a synchronously updated Boolean network is a function F = (f1, . . ., fn) : where each coordinate function fi describes how the future value of variable xi depends on the present values of all variables.All variables are updated at the same time (synchronously).Remark 2.3.The wiring diagram of any Boolean network is either strongly connected or it is composed of a collection of strongly connected components where connections between different components move in only one direction.

Dynamics of Boolean networks
Another directed graph associated with a BN is the state transition graph, also referred to as the state space.It describes all possible transition of the BN from one time step to another.The attractors of a BN are minimal sets of states from which there is no escape as the system evolves.An attractor with a single state is also called a steady state (or fixed point).In mathematical models of intracellular regulatory networks, the attractors of the model are often associated with the possible phenotypes of the cell.This idea can be traced back to Waddington [28] and Kauffman [29].For example, in a model of cancer cells, the steady states of the model correspond to proliferative, apoptotic, or growth-arrest phenotypes [30].Mathematically, a phenotype is associated with a group of attractors where a subset of the system's variables have the same states.These shared states are then used as biomarkers that indicate diverse hallmarks of the system.There are two ways to describe the dynamics of a Boolean network F : F n 2 → F n 2 , (i) as trajectories for all 2 n possible initial conditions, or (ii) as a directed graph with nodes in F n 2 .Although the first description is less compact, it will allow us to formalize the dynamics of coupled networks.Definition 2.5.A trajectory of a Boolean network F : We can see that T3 and T6 are periodic trajectories with period 2. Similarly, T1 and T4 are periodic with period 1.All other trajectories eventually reach one of these 4 states.
Definition 2.7.The state space of a (synchronously updated) Boolean network F : F n 2 → F n 2 is a directed graph with vertices in F n 2 and an edge from x to y if F (x) = y.
Example 2.8. Figure 2b shows the state space of the (synchronously updated) Boolean network from Example 2.4.
From the state space, one can easily obtain all periodic points, which form the attractors of the network.
Definition 2.9.The space of attractors for a Boolean network is the set 1.The subset D 1 (F ) ⊂ D(F ) of sets of exact size 1 consists of all steady states (also known as fixed points) of F .
2. The subset D r (F ) ⊂ D(F ) of sets of exact size r consists of all cycles of length r of F .
Equivalently, an attractor of length r is an ordered set with r elements, C = {c1, . . ., cr}, such that Remark 2.10.In the case of steady states, the attractor C = {c} may be denoted simply by c.
Example 2.11.The Boolean network from Example 2.4 has 2 steady states (i.e., attractors of length 1) and one cycle of length 2, which can be easily derived from its state space representation (Figure 2b).

Modules
In [25], a concept of modularity was introduced for Boolean networks.The decomposition into modules occurs on structural (wiring diagram) level but induces an analogous decomposition of the network dynamics, in the sense that one can recover the dynamics of the entire network from the dynamics of the modules.For this decomposition, a module of a BN is defined as a subnetwork that is itself a BN with external parameters in the subset of variables that specifies a strongly connected component (SCC) in the wiring diagram (see Example 2.12).More precisely, for a Boolean network F and subset of its variables S, we define the restriction of F to S to be the BN We note that f k i might contain inputs that are not part of S (e.g., when x k i is regulated by some variables that are not in S).Therefore, the BN F |S may contain external parameters (which are themselves fixed and do not possess an update rule).Given a F with wiring diagram W , let W1, . . ., Wm be the SCCs of W with pairwise disjoint sets of variables Si.The modules of F are then the restrictions to these sets of variables, F |S i .Further, the modular structure of F can be described by a directed acyclic graph Q = {(i, j) | Wi −→ Wj} by setting Wi −→ Wj whenever there exists a node from Wi to Wj (see Example 2.12).
Example 2.12.Consider the Boolean network with wiring diagram in Figure 3a.The restriction of this network to S1 = {x1, x2} is the 2-variable network ), which forms the first module (indicated by the amber box in Figure 3a), while the restriction of ), which forms the second module (indicated by the green module in Figure 3a).Note that the module F |S 1 , i.e., the restriction of F to S1, is simply the projection of F onto the variables S1 because W1 does not receive feedback from the other component.

Control via Modularity
In this section, we apply the modular decomposition theory described in the previous section and in [25] to make the control problem of Boolean networks more tractable.We show how the decomposition into modules can be used to obtain controls for each module, which can then be combined to obtain a control for the entire network.In this context, two types of control actions are generally considered: edge controls and node controls.For each type of control, one can consider deletions or constant expressions as defined below.The motivation for considering these control actions is that they represent the common interventions that can be implemented in practice.For instance, edge deletions can be achieved by the use of therapeutic drugs that target specific gene interactions, whereas node deletions represent the blocking of effects of products of genes associated to these nodes; see [31,32].
Once the modules have been identified, different methods for phenotype control (that is, control of the attractor space) can be used to identify controls in these networks.Some of these methods employ stable motifs [11], feedback vertex sets [13], as well as algebraic approaches [33,34,35].For our examples below, we will use the methods defined in [11,33,13] to find controls for the simple networks.
A Boolean network F = (f1, . . ., fn) : where U is a set that denotes all possible controls, defined below.The case of no control coincides with the original Boolean network, that is, F µ (x, 0) = F (x).Given a control µ ∈ U , the dynamics are given by x(t + 1) = F µ (x(t), µ).See [33] for additional details and examples of how to encode control edges and nodes in a Boolean network.
where ai is a constant in F2, encodes the control of the edge xi → xj, since for each possible value of µi,j ∈ F2 we have the following control settings: • If µi,j = 0, F µ j (x, 0) = fj(x1, . . ., xi, . . ., xn).That is, the control is not active.• If µi,j = 1, F µ j (x, 1) = fj(x1, . . ., xi = ai, . . ., xn).In this case, the control is active, and the action represents the removal of the edge xi → xj when ai = 0, and the constant expression of the edge if ai = 1.We use xi This definition can be easily extended for the control of many edges, so that we obtain , where e is the number of edges in the wiring diagram.Each coordinate, µi,j, of µ in F µ (x, u) encodes the control of an edge xi → xj.
encodes the control (knock-out or constant expression) of the node xi, since for each possible value of (µ we have the following control settings: This action represents the knock-out of the node xj.
This action represents the constant expression of the node xj.
. This action changes the Boolean function to its negative value.This case is usually not considered in the control search since it is biologically impractical to implement.
We note that the algebraic framework is versatile enough that we can encode any type of control, such as a combination of node and edge control at the same time.Definition 3.3.A control µ stabilizes a given Boolean network F : F n 2 → F n 2 at an attractor C ⊆ F n 2 when the resulting network after applying µ to F (denoted here as F µ ) has C as its only attractor.
For a Boolean network F , we let D(F ) denote the set of its attractors.Whenever a Boolean network F has more than one module we say that it is decomposable into its constituent modules F1, F2, • • • , Fm (m ≥ 2), and write Fm where the semi-product operation ⋊P i indicates the coupling of the subnetworks, as described in [25] (see Example 3.4 for an example of coupling).Furthermore, from the decomposition theory described in [25], the attractors of F are of the form C = C1 ⊕ C2 ⊕ • • • ⊕ Cn where Ci ∈ D(Fi) is an attractor of the subnetwork, for i = 1, . . ., n.The following theorem takes advantage of the modular structure of the network to find controls one module at a time.
Example 3.4.Consider the Boolean network in Example 2.12, with wiring diagram in Figure 3a.
which forms the first module (indicated by the amber box in Figure 3a) and let ) be the restriction of F to S2 = {x3, x4}, which forms the second module (indicated by the green box in Figure 3a).Here, e1 and e2 are external parameters of F |S 2 that take the place of x1 and x2, which is indicated by the coupling P = {x1 → e1, x2 → e2}.F is thus decomposable into F1 and F2 and we have F = F1 ⋊P F2.Theorem 3.5.Given a decomposable network F = F1 ⋊P F2, if µ1 is a control that stabilizes F1 in C1 (whether C1 is an existing attractor of F1 or a new one, after applying control µ1) and µ2 is a control that stabilizes F2 in C2 (whether C2 is an existing attractor of F2 or a new one, after applying control µ2), then µ = (µ1, µ2) is a control that stabilizes F in C = C1 ⊕ C2 provided that at least one of C1 or C2 is a steady state.
Proof.Let F µ 1 1 be the resulting network after applying the control µ1.Thus, the dynamics of be the resulting network after evaluating the external parameters of F2 in the states of C1 and applying the control µ2.Then, D(F C 1 ,µ 2 2 ) = C2.Thus, Therefore, For the last equality we used the fact that the product of a steady state and a cycle (or vice versa) will result in only one attractor for the combined network.The former is not always true in general because multiplying two attractors (of length greater than 1) might result in several attractors for the composed network due to the attractors starting at different states.It follows that there is only one attractor of F µ and that attractor is C1 ⊕ C2.Thus, F is stabilized by µ = (µ1, µ2) and we have D(F µ ) = C. Theorem 3.5 shows how the modular structure can be used to identify controls that stabilize the network in any desired state.In particular, we can use the modular structure of a network to find controls that stabilize a network at an existing attractor, which is often the case in biological control applications.We state this fact in the following corollary.
Corollary 3.6.Given a decomposable network F = F1 ⋊P F2, let C = C1 ⊕C2 be an attractor of F , where C1 ∈ D(F1) and C2 ∈ D(F C 1 2 ) and at least C1 or C2 is a steady state.If µ1 is a control that stabilizes F1 in C1 and µ2 is a control that stabilizes Remark 3.7.In Theorem 3.5, we required one of the stabilized attractors to be a steady state in order to be able to combine the controls from the modules.We can remove this requirement from Theorem 3.5 by using the following definition of stabilization for non-autonomous networks, which will guarantee that C1 and C2 can be combined in a unique way, resulting in a unique attractor of the whole network.Definition 3.8.A non-autonomous Boolean network is defined by where H : F m+n 2 → F n 2 and (g(t)) ∞ t=0 is a sequence with elements in F m 2 .We call this type of network non-autonomous because its dynamics will depend on g(t).We use H g to denote this non-autonomous network.
If H(g(t), y) = G(y) for some network G (that is, it does not depend on g(t)) for all t, then y(t+1) = H(g(t), y(t)) = G(y(t)) and this definition of attractors coincides with the classical definition of attractors for (autonomous) Boolean networks (Definition 2.9).Definition 3.9.Consider a controlled non-autonomous network given by y(t + 1) = F2(g(t), y(t), µ), where g(t) is a trajectory representation of an attractor C1 of an upstream network.We say that a control µ2 stabilizes this network, F C 1 2 (defined as in Definition 3.8), at an attractor C2 when the resulting network after applying µ2 (denoted here as ) has C2 as its unique attractor.For non-autonomous networks the definition of unique attractor requires that (g(t), y(t)) ∞ t=0 has a unique periodic trajectory up to shifting of t (which is automatically satisfied if C1 or C2 is a steady state). 2 ), which is not an original attractor of F C 1 2 .• Finally, the control µ = (µ1, µ2) stabilizes F at C = C1 ⊕ C2 = {0110}.Note that C is a new attractor for F .Theorem 3.5 uses the modular structure of a Boolean network to identify controls that stabilize the network in any desired attractor.In biological applications, the attractors typically correspond to distinct biological phenotypes (defined more rigorously in the next section) and a common question is how to force a network to always transition to only one of these phenotypes.For example, cancer biologists may use an appropriate Boolean network model with the two phenotypes proliferation and apoptosis to identify drug targets (i.e., edge or node controls), which force the system to always undergo apoptosis.The following example illustrates this specific control aspect, described in Corollary 3.6.• The edge control µ2 : x4 1 − → x3 (that is, the control that constantly expresses the edge from x4 to x3) stabilizes 2 ).Again, note that x3 1 − → x4 would be an alternative control.

Control via Modularity and Canalization
In addition to using the modular structure of the network, we can take advantage of the canalizing structure of the regulatory functions to identify control targets.We first review some concepts and definitions, and introduce the concept of canalization.
where ei is the ith unit vector.In that case, we also say f depends on xi.
The last line ensures that f actually depends on all n variables.
To account for partial canalization, we also define k-canalizing functions which were first introduced in [36].
where fC = fC (x σ(k+1) , . . ., x σ(n) ) is the core function, a Boolean function on n − k variables.When fC is not canalizing, then the integer k is the canalizing depth of f [37].Note that an n-canalizing function (i.e., a function where all variables become eventually canalizing) is also called a nested canalizing function (NCF).
When f is not canalizing (i.e., when k = 0), we simply have pC = f .Definition 4.6.Given a Boolean function f (x1, . . ., xn) represented as in Equation 3, we call the extended monomials Mi the layers of f and define, as in [38], the layer structure as the vector (k1, . . ., kr), which describes the number of variables in each layer.Note that f is nested canalizing if and only if k1 Remark 4.7.Here we note the following important properties of layers of canalization.
(a) Theorem 4.5 shows that any Boolean function has a unique extended monomial form given by Equation 3, in which the variables are partitioned into different layers based on their dominance.Any variable that is canalizing (independent of the values of other variables) is in the first layer.Any variable that "becomes" canalizing when excluding all variables from the first layer is in the second layer, etc.All remaining variables that never become canalizing are part of the core polynomial.The number of variables that eventually become canalizing is the canalizing depth of the function.NCFs are exactly those functions where all variables become eventually canalizing (note not all variables of an NCF must be in the first layer).
(b) While variables in the same layer may have different canalizing input values, they all share the same canalized output value, i.e., they all canalize a function to the same output.On the other hand, the outputs of two consecutive layers are distinct.Therefore, the number of layers of a k-canalizing function, expressed as in Definition 4.4, is simply one plus the number of changes in the vector of canalized outputs, (b1, b2, . . ., b k ).• x1 canalizes f to 0 if it receives its canalizing input 0.
• If this is not the case, x2 = 0 canalizes f to 1.
While finding the layer structure of a Boolean function is an NP-hard problem, there exist several algorithmic implementations [39].
In meaningful biological networks, the attractors correspond to phenotypes.Typically a small subset of all Boolean variables is used to define phenotypes.Definition 4.9.Given a Boolean network F with attractor space D(F ) and phenotype-defining variables P ⊂ {x1, . . ., xn}, we associate the same phenotype to all attractors C ∈ D(F ) that are identical in P. The states in P will be called markers of the phenotype.
Suppose F = F1 ⋊P F2 is a decomposable network, and that there is a phenotype that depends on variables in F2 only (that is, all markers of the phenotype are part of F2), and that we wish to control the phenotype through F2.The most straightforward approach is to set the variables that the phenotype depends on to the appropriate values that result in the desired phenotype.However, such intervention may not be experimentally possible.Instead, we can exploit the canalizing properties of the functions corresponding to the nodes connecting the modules F1 and F2 to identify control targets.
By Theorem 4.5, the variables of any Boolean update function can be ordered by importance/dominance, based on which layer they appear in (with variables in outer layers being more important).Thus, once we control a variable in a certain layer (by setting it to its canalizing input value), any further control of variables in less dominant layers will have no effect on the function (and thus on the network).We state this fact in the following lemma.Lemma 4.10.Suppose F = F1 ⋊P F2 is a decomposable network.Suppose further that only one node x ∈ F2 with update function fx is regulated by nodes in F1.If fx is canalizing with r layers, let ℓ ∈ {1, . . ., r} be the lowest (i.e., most important) layer of fx, which contains nodes from F1.If all regulators of x from F1 appear in the core polynomial, we set ℓ = r + 1.Then, setting y ̸ ∈ F1 to its canalizing value decouples the systems F1 and F2, as long as y appears in a layer < ℓ.When a node in module F m is regulated by a node in module F k (indicated by the red edge in the directed acyclic graph), the phenotype may no longer depend on nodes in F m , in order for module F i to be removable from the control search.
Proof.The lemma is a direct consequence of Theorem 4.5.If y receives its canalizing input and is in a more important layer of fx than all variables in F1, then none of these variables can affect fx anymore.Thus, controlling y to receive its canalizing input eliminates the link between F1 and F2.
The modularity approach in [25] yields an acyclic directed graph after one collapses each module into a single node.This endows a natural partial ordering on the collection modules of a network where one module precedes the next if there is a path for every node in the first module which ends in the second module.While not all modules are comparable, we can speak of chains of modules which consist of subsets of the partial ordering which are totally ordered.Furthermore, one can rank the modules based on the percentile scores (i.e., rank module k out of m modules).This type of ranking has been studied in [40], where it was shown that the importance of the modules is strongly correlated with the aggressiveness of mutations occurring within those modules and the effectiveness of interventions.
Theorem 4.11.Suppose F = F1 ⋊P 1 F2 ⋊P 2 • • • ⋊P n−1 Fm is a decomposable network.If for some i < j, (i) only one node x ∈ Fj with update function fx is regulated by nodes in Fi, and (ii) fx is a canalizing function, which possesses none of the variables from Fi in its most important layer, and (iii) the phenotype of interest only depends on variables in Fj as well as modules F k , for which any chain containing Fi and F k also contains Fj. (see Fig. 4 for an example) then the module Fi can be excluded from the control search by setting any node y / ∈ Fi to its canalizing input, as long as this node appears in a more dominant layer of fx than all inputs of fx that are part of Fi.
Proof.By Lemma 4.10, setting y to its canalizing value results in decoupling Fi and Fj.Fi will no longer have any effect on Fj, and thus, due to condition (iii), on the phenotype of interest.Fi can therefore be removed from the control search.Theorem 4.11 is illustrated in Fig. 5.Note that node y can be in Fj or some other module as in the figure.In this theorem, we assumed that none of the variables of Fi are in the most dominant layer in the update rules of variables in Fj.If some of the variables of Fi are in the most dominant layer, we can still remove module Fi from the control search using an edge control, as shown in the following theorem.Theorem 4.12.Suppose F = F1 ⋊P 1 F2 ⋊P 2 • • • ⋊P n−1 Fm is a decomposable network.If for some i < j, (i) only one node x ∈ Fj with update function fx is regulated by nodes in Fi, and (ii) fx is a canalizing function with some variables from Fi in its most important layer, and

Phenotype
Figure 5: Control via modularity and canalization.Once the network is decomposed into modules F 1 , . . ., F m , we can override the effect of module F i by the using another module (F k in this case) whose variables are inputs of f x that are located in a higher canalizing hierarchy than the layers containing the variables of F i .
(iii) the phenotype of interest only depends on variables in Fj as well as modules F k , for which any chain containing Fi and F k also contains Fj. (see Fig. 4 for an example) then the module Fi can be excluded from the control search by applying an edge control to any input in the most dominant layer of fx.
Proof.Let y ∈ Fi such that y ∈ supp (fx), and that y is located in the most dominant layer fx.Then, setting y to its canalizing value results in decoupling the subnetworks Fi and Fj.Thus, Fi will no longer have any effect on Fj and thus it can be removed from the control search.
Remark 4.13.The method can be extended to the case when Fi and Fj are connected via multiple nodes (that is, condition (i) in Theorem 4.11 and Theorem 4.12 can be relaxed.In that case decoupling is achieved through the same procedure presented above, applied to each node in Fj with regulators from Fi.
We further note that condition (ii) in the theorems above is generally very restrictive as only a small proportion of Boolean functions in n > 3 variables is canalizing or even nested canalizing.However, as shown in [26], most biological Boolean network models are almost entirely governed by nested canalizing functions.
To showcase these methods, we will now decompose a published Boolean network model into its modules, and then identify the minimal set of controls for the entire network by exploiting the canalizing structure of the regulatory functions within the modules.The identified set of controls will force the entire system into a desired attractor.
Example 4.14.We consider a Boolean network model for the blood cancer large granular lymphocyte (T-LGL) leukemia, which was published in [27].T-LGL leukemia is a clonal hematological disorder characterized by persistent increases of large granular lymphocytes in the absence of reactive cause [41].The wiring diagram of this model is depicted in Figure 6a.This network has 16 nodes and three nontrivial modules (highlighted by the amber, green, and gray boxes in Figure 6a).The control objective here is to identify control targets that lead the system to programmed cell death.In other words, we aim to direct the system into an attractor that has the marker apoptosis ON.
Since the model has three nontrivial modules, the approach described in Section 3 would require us to identify control targets for three modules.However, an exploitation of the canalizing structure and common sense reveals that we do not need to control every module to ensure apoptosis, the desired control objective.First, irrespective of canalization, the module highlighted in gray in Figure 6a does not affect the phenotype apoptosis.Therefore, we can focus on the modules "upstream" of apoptosis (i.e., the amber and green modules in Figure 6a).In this case, we will apply Theorem 4.12 to identify control targets for this model.We note that the edges from the upstream module (amber box in Figure 6a) to the downstream module (green box in Figure 6a) all end in the node DISC.Therefore, we will investigate the canalizing properties of the regulatory function of DISC (see Figure 6b), fDISC = Ceramide ∨ (F as ∧ F LIP ).
Using the approach described in [39], we find that fDISC has two canalizing layers, L1 = {Ceramide} and L2 = {F as, F LIP } and its layering structure is given by (see Figure 6c) fDISC = (Ceramide + 1)[(F as)(F LIP + 1) + 1]+1 We note that the only variable in the most important canalizing layer, Ceramide, is in the upstream module.Thus, we can decouple the modules via an edge control on the connection between the upstream and downstream modules.That is, the constant expression of the edge from Ceramide to DISC will decouple the two modules and will lead to constant expression of DISC.We can check that this control is effective at stabilizing the system in the desired attractor and the control set obviously has minimal size.
In summary, in this example we used an edge control to decouple the upstream and downstream modules and then identified a control target in the downstream module which contains the markers of the phenotype of interest.

Conclusion
Model-based control is a mainstay of industrial engineering, and there is a well-developed mathematical theory of optimal control that can be applied to models consisting of systems of ordinary differential equations.While this model type is also commonly used in biology, for instance in biochemical network modeling or epidemiology and ecology, there are many biological systems that are more suitably modeled in other ways.Boolean network models provide a way to encode regulatory rules in networks that can be used to capture qualitative properties of biological networks, when it is unfeasible or unnecessary to determine kinetic information.While they are intuitive to build, they have the drawback that there is very little mathematical theory available that can be used for model analysis, beyond simulation approaches.And for large networks, simulation quickly becomes ineffective.
The results in this paper, building on those in [25], can be considered as a contribution to a mathematical control theory for Boolean networks, incorporating key features of biological networks.There are many open problems that remain, and we hope that this work will inspire additional developments.
Our concrete contributions here are as follows.The modularization method makes the control search far more efficient and allows us to combine controls at the module level obtained with different control methods.For example, methods based on computational algebra [33,35] can identify controllers that can create new (desired) steady states, which other methods cannot.Feedback vertex set [14,13] is a structure-based method that identifies a subset of nodes whose removal makes the graph acyclic.Stable motifs [11] are based on identifying strongly connected subgraphs in the extended graph representation of the Boolean network.Other control methods include [42,43,44].We can use any combination of these methods to identify the controls in each module.

Definition 2 . 1 .Definition 2 . 2 .
The wiring diagram of a Boolean network F = (f1, . . ., fn) :F n 2 → F n2 is the directed graph with vertices x1, . . ., xn and an edge from xi to xj if fj depends on xi.That is, if there exists x ∈ F n 2 such that fj(x) ̸ = fj(x + e i ), where e i is the ith unit vector.The wiring diagram of a Boolean network is strongly connected if every pair of nodes is connected by a directed path.That is, for each pair of nodes xi, xj in the wiring diagram with xi ̸ = xj there exists a directed path from xi to xj (and vice versa).In particular, a one-node wiring diagram is strongly connected by definition.

Example 2 . 6 .Figure 2 :
Figure 2: Wiring diagram and state space of the Boolean network in Example 2.4-2.11.(a) The wiring diagram encodes the dependency between variables.(b) The state space is a directed graph with edges between all states and their images.This graph therefore encodes all possible trajectories.

Figure 3 :
Figure 3: Boolean network decomposition into modules.(a) Wiring diagram of a non-strongly connected Boolean network where the non-trivial modules are highlighted by amber and green boxes.(b) Directed acyclic graph describing the corresponding connections between the nontrivial modules.

Definition 3 . 2 (
Node Control).Consider the node xi in the wiring diagram W .The function

Definition 4 . 1 .
A Boolean function f (x1, . . ., xn) is essential in the variable xi if there exists an

k i j=1 (
xi j +ai j ) is a nonconstant extended monomial, pC is the core polynomial of f , and k = r i=1 ki is the canalizing depth.Each xi appears in exactly one of {M1, . . ., Mr, pC }, and the only restrictions are the following "exceptional cases":

Figure 4 :
Figure 4: Example of a modular directed acylic graph structure to illustrate Condition (iii) in Theorems 4.11 and 4.12.(a) Module F i can be removed from the control search as long as conditions (i) and (ii) are satisfied, and the phenotype of interest depends only on any subset of variables that are part of the blue modules.(b)When a node in module F m is regulated by a node in module F k (indicated by the red edge in the directed acyclic graph), the phenotype may no longer depend on nodes in F m , in order for module F i to be removable from the control search.

Figure 6 :
Figure 6: (a) Wiring diagram of the T-LGL model, published in [27], which describes the mechanisms that regulate apoptosis.The non-trivial modules are highlighted by amber, green, and gray boxes.(b) The regulatory inputs of the node DISC.(c) Writing the regulatory function corresponding to node DISC in its standard monomial form (Theorem 4.5) reveals its canalizing structure.